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ABSTRACT 

The artificial viscosity is reconsidered in smoothed particle hydrodynamics to 
prevent inter-particle penetration, unwanted heating, and unphysical solutions. 
The coefficients in the Monaghan's standard artificial viscosity are considered 
as time variable, and a restriction on them is proposed such that avoiding the 
undesired effects in the subsonic regions. The shock formation in adiabatic and 
isothermal cases are used to study the ability of this modified artificial viscosity 
recipe. The computer experiments show that the proposal appears to work and 
the accuracy of this restriction is acceptable. 

Subject headings: Hydrodynamics - methods: numerical - ISM: evolution 

1. Introduction 

Gas dynamical processes are believed to play an important role in the evolution of astro- 
physical systems on all length scales. A ubiquitous process in astrophysical fluid dynamics 
is the behavior of gases subjected by shock waves. These are almost the normal, rather than 
exceptional, type of astrophysical fluid flows that happen very often in cases of astrophys- 
ical interest. A number of examples are the onset of spiral structures and shock fronts in 
accretion disks around compact objects (Lanzafame et al. 2006), molecular cloud formation 
in shock-compressed layers (Koyama and Inutsuka 2000), the rapid and strong increase in 
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pressure that can be produced by interaction of the fast wind and the slow envelope of a 
planetary nebulae (Gurzadyan 1997), and so on. 

A powerful gridless particle method that was invented to solve complex fluid-dynamical 
problems in astrophysics is smoothed particle hydrodynamics (SPH) (Lucy 1977, Gingold 
and Monaghan 1977). The SPH method has a number of attractive features such as its 
low numerical diffusion in comparison to grid based methods. Whilst the SPH originally 
developed for compressible flows, it has been extended to deal with free-surface incompress- 
ible fluids (e.g., Monaghan 1994, Ellero et al. 2007). A recent worthy review of the SPH 
methodology and its applications can be found in Monaghan (2005). 

An adequate scenario for SPH application is the unbound astrophysical problems, espe- 
cially the behavior of gases subjected to compression (shock waves). The magnitude of the 
viscous term dose not affect the net shock-jump conditions, thus, many numerical schemes 
implicitly or explicitly incorporate the trick of artificial viscosity (AV) for halting the ever- 
growing steepening tendency produced by nonlinear effects. The first use of viscosity in SPH 
equations was by Lucy (1977) who introduced an AV to present a slow build-up of acoustic 
energy from integration errors in an SPH simulation. The AV is designed to allow shock 
phenomena to be simulated, or simply to stabilize a numerical algorithm. Indeed two forms 
of AV are applied in SPH equations, the bulk and the von Neumann-Richtmeyer viscosity, 
respectively. They prevent inter-particle penetration, allow shocks to form and damp the 
post-shock oscillations. Although many different functional forms for the AV have been pro- 
posed (see, e.g., Liu & Liu 2003), all contain problem sensitive parameters that are often 
set in a somewhat arbitrary manner. Here we try firstly to provide a guide to the relevant 
prescriptions for some AV that have been considered before. 

A more effective viscosity which may conserves linear and angular momentum was sug- 
gested by Monaghan and Gingold (1983). They devised a viscosity by simple arguments 
about its form and its relation to gas viscosity. The viscous term between two particles, de- 
noted by n a b, is added to the pressure term of SPH equations. When two particles approach 
each other, the AV produces a repulsive force between them. When they recede from each 
other the force is attractive (Monaghan 1989). 

An undesirable aspect of the Monaghan's standard AV is that it can introduce con- 
siderable shear viscosity into the flow. This problem may be reduced by a bulk viscosity 
in a general-purpose code of Hernquist and Katz (1989) for evolving three-dimensional, 
self-gravitating fluids in astrophysics. Clearly the artificial viscous dissipation increases 
the Reynolds number of a flow, artificially, with the result that, for example, the Kelvin- 
Helmholtz shear instabilities are heavily diffused. In this way a von Neumann stability 
analysis of the SPH equations along with a critical discussion of various parts of the al- 
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gorithm was investigated by Balsara (1995). He suggested reducing viscous dissipation by 
multiplying U ab by a symmetric factor. 

Another problem of viscosity is that although AV is approximately successful for han- 
dling shocks but it can be too large in other parts of the flow. For this purpose, a switch 
to reduce the AV away from shocks was given by Morris and Monaghan (1997, hereafter 
MM97). They introduce the idea of time- varying switch which fits more naturally with a 
particle formulation. Each particle has a viscosity parameter which evolves according to a 
simple source and decay equation. The source causes the switch to grow when the particle 
enters a shock and the decay term causes it to decay to a small value beyond the shock. 

In the present study we combine the Monaghan's standard AV with the time-varying 
coefficients like that of MM97's switch, we propose a modification to the time-dependent 
AV prescription designed by MM97 so that the aim is to maintain the shock-simulating 
capabilities of SPH. For this purpose we optimize the source term so that the AV is restricted 
to supersonic velocities in regions that are under compression. Section 2 is devoted to the 
SPH method and the suggested changes in the AV. These modifications have been tested in 
§3 for the one dimensional adiabatic and isothermal shock problems. Finally, a summary 
with conclusion is given in section 4. 



2. The SPH methodology 

The SPH was invented to simulate nonaxisymmetric phenomena in astrophysics (Lucy 
1977, Gingold & Monaghan 1977). In this method, fluid is represented by iV discrete but 
extended/smoothed particles (i.e. Lagrangian sample points). The particles are overlapping, 
so that all the involved physical quantities can be treated as continuous functions both in 
space and time. Overlapping is represented by the kernel function, W a b = W(r a — r^, h a b), 
where h a b = (h a +hb)/2 is the mean smoothing length of two particles a and b. The continuum 
equations are (Monaghan 1992) 

Pa = ^2m b W ab (1) 

b 

^ = - E m ^~2 + + U ^ V « W « b ( 2 ) 

d -ir = \Y. m ^-2 +~2+ n ab )v afc • v a w ab (3) 

at b P a P h 

where v a ;, = v a — v;,, the AV between particles a and b is presented by U a b, and the other 
notations have their usual meanings. 
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2.1. Viscosity 

In order to simulate the shock problems of hydrodynamics, special treatments or meth- 
ods are required to allow the algorithms to be capable of modeling shock waves, or else the 
simulation will develop unphysical oscillations in the numerical results around the shocked 
regions. A shock wave is not a true physical discontinuity, but a very narrow transition zone 
whose thickness is usually in the order of a few molecular mean free paths. Application of 
the conservation of mass, momentum, and energy conditions across a shock front requires 
the simulation of transformation of kinetic energy into heat energy. Physically this energy 
transformation can be represented as a form of viscous dissipation. This idea leads to the 
development of the von Neumann-Richtmyer AV that is needed only to be present during 
material compression, 

= U(Ax) 2 p(V.v) 2 V.v<0 

1 \0 V.v > 0, 1 J 

where a± is an adjustable non-dimensional constant. It is found that adding the following 
linear AV term, 

= fa 2 (Ax)cp(V.v) V.v<0 

2 \0 V.v > 0, 1 ' 

where c is the sound speed and a 2 is another non-dimensional constant (Liu and Liu 2003), 
has the advantage of further smoothing the oscillations that are not totally dampened by 
the quadratic AV term, equation (J3J. 

The most widely used AV in SPH is that of Monaghan (1989), which not only provides 
the necessary dissipation to convert kinetic energy into heat at the shock front, but also pre- 
vent unphysical penetration for particles approaching each other. The detailed formulation 
is as follows: 

n ab = \ ^ > 11 V >*- T <* < u ' (6) 

1 0, otherwise, 



where p ao = |(p a + pb) is an average density, a ~ 1 and (3 ~ 2 are the artificial coefficients, 
and p a b is defined as its usual form 

p — Vab ' rab - (7) 

Kb r 2 Jh 2 ab + T] 2 

with 77 ~ 0.1 and h a b = \{h a + h ). The signal velocity, v S i 9 , is 

v sig = ^(c a + c b ) (8) 



where c a and % are the sound speed of particles. 
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Since the Monaghan type AV introduces a shear viscosity into the flows especially in 
regions away from the shock, an AV depending on the divergence of the velocity field was 
employed by Hernquist and Katz (1989), 

IU = % + % (9) 

pl pI 



q a = J ah a p a c a \ V.v| + (3h 2 a p a \V .v\ 2 a V.v < ^ 



where 

Qn = < 

otherwise. 
This AV can lead to spurious results if particles a and b are receding; i.e., such terms will 
contribute viscous artificial cooling. To avoid this problem, Hernquist and Katz (1989) set 
Ii a b = if v a fc.r a b > 0; a decision which may, in some cases, add small amounts of shear 
viscosity. In this approach, determining Ii ab by referring to binary interaction, effectively, 
gives an idea about the molecular bulk fluid descriptions. 

Other view on the Monaghan's standard AV was also proposed by Balsara (1995) who 
developed a modification of the AV term that approximately sets the dissipation to zero in 
regions of pure shear, while leaving it unaffected in regions of compression. This is done by 
suggestion 

rU = ( % + ^ J {-ap ab + (3p 2 ab ), (11) 

\Pa Pb J 

where a switch in the form of a multiplicative factor (p ab — > p a bfab/c a b) was introduced. The 
factor f a b is the average of f a and fb, where 





V.v 


1 o, 


|V.v: 




V x v| 


a + lO- 4 C a /h a 



The function f a acts as a switch, approaching unity in regions of compression (|V.v| a >> 
|V x v| a ) and vanishing in regions of large vorticity (|V x v| >> |V.v| ). Consequently, 
this AV has the advantage that it is suppressed in shear layers. 

Meanwhile, MM97 carried the process further by introducing the idea of time-varying 
coefficients which fits more naturally with a particle formulation. In this method, the AV is 
given by 

K a bVgigV a b " Y a b . . 

Waft = — -; — ; , [i-O) 

Pab\Yab\ 

where K ab = (K a + K b )/2 is a variable switch which should change with time according to 
the conditions the particle is in, becoming large at shocks, but relaxing back to a small value 
when the flow is calmer. Each particle a has a viscosity switch which evolves according to a 
simple source, 

S a = max(-V-v a ,0), (14) 
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and decay equation, 

dK K — K ■ 

-7- = + So, 15 

at r a 

such that in the absence of source term <S a , the switch K a decays to a value K min over a time 
scale r a = h a /(Cv sig ). The dimensionless parameter C has a value 0.1 < C < 0.2. The source 
causes the switch to grow when the particle enters a shock and the decay term causes it to 
decay to a small value beyond the shock. 



2.2. Restricted AV 

In SPH the formation of shocks is mostly an effect of the von Neumann-Richtmeyer 
viscosity. To avoid the undesirable effects in the subsonic region, we propose to restrict the 
use of it to those regions. This modification will allow shocks to form, and prevent inter- 
particle penetration at supersonic velocities. When the gas has reached equilibrium, the 
pressure force prevents further compression. This process is not possible with source term 
( rHj) because of supersonic relative velocities of the particles at different regions of the shock. 
The AV as outlined by equation (fl3|) will decelerate and heat the gas, and prevent enough 
compression. To weaken this problem, we use the following lemma. 

Since p a oc h~ u where v is the dimension of the problem, the time derivative for any 
particle can be written as 

Pa = 7 • (16) 

"a 

This relation can be used to decide whether a particle follows the fluid or if the source term 
of AV is necessary. If two particles are approaching each other at velocity exceeding the 
signal velocity, that is if 

r ab (— - ^) > v sig , (17) 

Pa i^a 

AV is necessary to prevent inter-particle penetration. We therefore propose to restrict the 
source term of AV with defined frequency as 

u; a = ^-^-^. (18) 
Pa K r ab 

A good way to do this for a particle is to use the restricted source term S a , instead of 
equation (TH1) . as follows 

S a = max(-V.v a ,^ a ,0). (19) 
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Base on a method similar MM97, we use a and (3 in the standard AV, equation ([6]), in 
form of variables with respect to time, 



da a 
dt 

and 



h (20) 



= — + (21) 

respectively, where the parameters z a and are chosen to regulate the effect of source term 
such that the heat production and post-shock oscillations are controlled in the numerical 
simulations. 



3. Shock formation test 

3.1. Analytical formulation 

As outlined before, an extremely important problem is the behavior of gases subjected 
to compression waves. This happens very often in the cases of astrophysical interests. For 
example, a small region of gas suddenly heated by the liberation of energy will expand into 
its surroundings. The surroundings will be pushed and compressed, thus, a shock front is 
formed. Conservation of mass, momentum, and energy across a shock front is given by the 
Rankine-Hugoniot conditions (Dyson and Williams 1997) 

Pivi = P2V2 (22) 

p x v\ + KpJ = p 2 v 2 2 + Kpl (23) 

\v\ + -^-KpJ- 1 = \v\ + ^-rKpV + Q (24) 
I 7— 1 I 7— 1 

where the equation of state, p = Kp 1 , is used. In adiabatic case, we have Q = 0, and for 
isothermal shocks, we will set 7 = 1. 

We would interested to consider the collision of two gas sheets with velocities vq in the 
rest frame of the laboratory. In this reference frame, the post-shock will be at rest and the 
pre-shock velocity is given by v\ = Vq + 1)2, where t>2 is the shock front velocity. Combining 
equations (!22l) - (l24l) . we have 

r b I b~ 2 " w M n 2 " , s 

V2 = a [-- + y 1 + - + (7 - 1)(^ - q)] (25) 
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where a = jKpJ 1 is the sound speed, M = v /a is the Mach number, b and q are defined 

as 

6= — Mo+ ^ 9= ^' (26) 

respectively. Substituting (I25I) into equation (1221 . density of the post-shock is given by 

p 2 =Pi{l + M ° (27) 



[-! + V 1 + T + (7-l)(T-?)] 



3.2. Simulation results 

Considering the head-on collision of two gas sheets, the particles with a positive x- 
coordinate are given an initial negative velocity of Mach 5, and those with a negative re- 
coordinate are given a Mach 5 velocity in the opposite direction. The chosen physical scales 
for length and time are [/] = 3.0 x 10 16 m and [t] = 3.0 x 10 13 s, respectively, so the velocity 
unit is approximately 1 km.s -1 . The gravity constant is set G — 1 for which the calculated 
mass unit is [m] = 4.5 x 10 32 kg. There is considered two equal one dimensional molecular 
sheets with extension x = 0.1 [Z], which have initial uniform density and temperature of 
~ 4.5 x 10 9 m -3 and ~ 10K, respectively. 

In adiabatic shock, with M = 5, the post-shock density must be 2.9, which is obtained 
from analytic solution (1271) with Q = and 7 = 2. On the other hand, in isothermal shock 
with Mq = 5, the post-shock density must be 26.9 (i.e., from Eq. [27] with 7 = 1). The 
simulation results of adiabatic and isothermal shocks which use the MM97's AV, equation 
f fT3|) . are shown in Fig. 1. The post-shock density has oscillations for both cases of the 
adiabatic and isothermal shocks. Since the thermal energy in adiabatic shock is restrained 
in the system, the oscillations grow further in it as shown in Fig. 1. Adding the linear 
AV term ([5]) can smooth the oscillations, thus, we use the parameter z a in the a viscosity 
equation (1201) . The best values of the z a for adiabatic and isothermal shocks are 2.0 and 3.5, 
respectively. The value of z a in adiabatic case is small for the reason that the oscillations 
are greater than isothermal case. 

As shown in Fig. 1, the post-shock density of adiabatic case is approximately equal to 
the analytical value, while the difference in isothermal case is more distinguishable. This 
is because of extremely increasing of the pressure in the compression region that is greatly 
produced by the /3-viscous term of AV in equation To remove this problem we propose 
using of the parameter zp in the (3 viscosity equation fT2T]) . Choosing the values 0.5 and 0.01 
for zp in adiabatic and isothermal case, respectively, can regulate the post-shock density as 
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shown in Fig. 2. These relative values are because the pressure produced in the post-shock 
region of isothermal case is artificially greater than the same pressure at adiabatic case. 

4. Summary and conclusions 

In SPH method two forms of AV is necessary; the von Neumann- Richtmeyer viscosity 
and linear viscous term. In this work we revisited the relevant prescriptions for some AV, 
which was proposed to consider these two forms of viscosity in SPH methodology. As a good 
idea, MM97 introduced the time-varying switch which fits more naturally with the particle 
formulation. The shock formations in adiabatic and isothermal cases have been performed 
to test the effect of different AVs. 

The simulation results in which the MM97's AV was used, show that the post-shock 
density has oscillations for both cases of the adiabatic and isothermal shocks. The oscillations 
in adiabatic shock grow further because the thermal energy in this case is restrained in the 
system. On the other hand, the post-shock density of adiabatic case is approximately equal 
to the analytical value, while the difference in isothermal case is more distinguishable. 

To partly overcome some of these problems, we proposed a modification of the AV. 
Based on the method of MM97, we used time-variable a and f3 in the standard AV. In this 
approach, the AV is used only in supersonic regions where the gas is not under compression, 
thus, the source term is restricted. Using this source multiplied by two regulating parameters 
(z a and zp) can virtually eliminate some problematic effects of using an artificial dissipation 
in SPH capability of shock simulation. Computer experiments have been performed to test 
the abilities of the proposal, and the results show that the accuracy of these restrictions on 
AV is acceptable. Although there are some satisfying results from shock formation tests, the 
proposed AV should be communicated to the astrophysical community, so that others can 
try it. There must be other tests to complete the idea. 
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Fig. 1. — The density of (a) adiabatic and (b) isothermal shocks in which the MM97's AV, 
equation (fl3l) . was used. 




Fig. 2. — The density of (a) adiabatic and (b) isothermal shocks in which the time variable 
a and (3 as outlined by equations (|20|) and (l2Tj) . were used. 



